%% Explore some data for cardiac output
% clear all; 
close all; clc;

DOSAVE = 0;
dirout = '../outputs/LVEFanalysis/';

% Load ptstruct
load '../data/ptstruct.mat';
N = length(ptstruct);

%% determine how many patients have more than one ICU admission
nicuadm = zeros(N,1);
npt_pulsep = 0;
npt_hr = 0;
npt_pulsep_hr = 0;
n = 0;
for i = 1 : N 
    nicuadm(i) = ptstruct(i).nicus(1);
    if nicuadm(i) == 1
        if ~isempty(ptstruct(i).data{1}.pulsep)
            npt_pulsep = npt_pulsep + 1;
        end
        if ~isempty(ptstruct(i).data{1}.hr)
            npt_hr = npt_hr + 1;
        end
        if ~isempty(ptstruct(i).data{1}.pulsep) && ~isempty(ptstruct(i).data{1}.hr)
            npt_pulsep_hr = npt_pulsep_hr + 1;
        end
    end
end

% represent results
[yicuadm,xicuadm] = hist(nicuadm,1:10);
fig(1) = figure;
% build axes positions
hBig = [subplot(121) subplot(122)];         %# create subplots
posBig = get(hBig, 'Position');             %# record their positions
delete(hBig)                                %# delete them
posSmall{1} = [0.275 0.63 0.16 0.24];
posSmall{2} = [0.717 0.63 0.16 0.24];
% create axes (big/small)
hAxB(1) = axes('Position',posBig{1});
hAxB(2) = axes('Position',posBig{2});
hAxS(1) = axes('Position',posSmall{1});
hAxS(2) = axes('Position',posSmall{2});
% plot
bar(hAxB(1), xicuadm, yicuadm, 'b');
bar(hAxB(2), 1:3, [npt_hr npt_pulsep npt_pulsep_hr], 'b');
bar(hAxS(1), xicuadm, yicuadm/N, 'r');
bar(hAxS(2), 1:3, [npt_hr npt_pulsep npt_pulsep_hr]/yicuadm(1), 'r');
set(hAxB(2),'XTickLabel',{'HR','PP','HR and PP'});
set(hAxS(2),'XTickLabel',{'HR','PP','HR and PP'});
%set(gca,'Xlabel','Nr ICU admissions','Ylabel','Count (patients)');
%# set axes properties
box off
set(hAxS , 'Color','none', 'XAxisLocation','bottom', 'YAxisLocation','left', 'box','off');
set(hAxS(1) , 'xlim',[.5 5.5]);
xlabel(hAxB(1),'Nr ICU admissions'); ylabel(hAxB(1),'Count (patients)');
xlabel(hAxB(2),'Patients with data'); ylabel(hAxB(2),'Count (patients)');
ylabel(hAxS(1),'Freq.');
ylabel(hAxS(2),'Freq.');
if DOSAVE
    fileout = sprintf('%snicuadm.png',dirout);
    print(fig(1),'-dpng',fileout);
end
%% study LVEF Data

pt_ind = [];
npt = 0;
los = [];
lvef_class{6} = [];
for i = 1 : N 
    nicuadm(i) = ptstruct(i).nicus(1);
    if ptstruct(i).nicus(1) == 1
        if ~isempty(ptstruct(i).data{1}.hr)
            npt = npt + 1;
            % pt_ind = [pt_ind i];
            los = [los ptstruct(i).los];
            clvef = ptstruct(i).clinics.lvef; % lvef category
            if (clvef > 4 || clvef < 0)
                lvef_class{6} = [lvef_class{6} i];
            else %if ptstruct(i).clinics.vp == 0 % include vasopressors?
                lvef_class{clvef+1} = [lvef_class{clvef+1} i];
            end     
        end
    end
end

fig(2) = figure; subplot(121); hist(los,1:2:100); xlabel('LoS (days)'); ylabel('Count (patients)');
xlim([0 100])
subplot(122); h = cdfplot(los); xlabel('LoS (days)'); xdata = get(h,'XData'); ydata = get(h,'YData');
xlim([0 100]); title(''); grid off
% plot percentiles
hold on;
pct = [95 98 99];
for p = pct
    [~, ip] = min(abs(xdata-prctile(los,p)));
    plot([0 xdata(ip) xdata(ip)], [ydata(ip) ydata(ip) 0], 'r');
end
if DOSAVE
    fileout = sprintf('%slos.png',dirout);
    print(fig(2),'-dpng',fileout);
end

% % Confirm that we assigned all patients
npt1 = 0;
nlvef = zeros(5,1);
for c = 1 : 5
    npt1 = npt1 + length(lvef_class{c});
    nlvef(c) = length(lvef_class{c});
end

fig(3) = figure;
bar(1:5,nlvef); set(gca,'XTickLabel',{'0','1','2','3','4'});
xlabel('LVEF'); ylabel('Count (patients)');
title(sprintf('Number of patients in each LVEF category (N = %d)',sum(nlvef)));
if DOSAVE
    fileout = sprintf('%slvef_classes.png',dirout);
    print(fig(3),'-dpng',fileout);
end

%% Analyse mean pulse pressure for different patient categories

test_time = [1 2 3 6 12 24 48 72 Inf];

for j = test_time
    lvef_class_pp = [];
    ISMEAN = 1;
    Npt = zeros(5,1);
    for c = 1: 5
        for i = 1 : length(lvef_class{c})
            ind = find(ptstruct(lvef_class{c}(i)).data{1}.hr(:,1) < j*60);
            if ~isempty(ind)
                vec = ptstruct(lvef_class{c}(i)).data{1}.hr(ind,2);
                if ISMEAN
                    vec = mean(vec(vec>0));
                else
                    vec = vec(vec>0);
                end
                Npt(c) = 1 + Npt(c);
                lvef_class_pp = [lvef_class_pp; vec (c-1)*ones(size(vec))];
            end
        end
    end
    classLVEF = 0:4;
    ppmean = zeros(1,5);
    
    fig(4) = figure(101);
    boxplot(lvef_class_pp(:,1),lvef_class_pp(:,2));
    for i = classLVEF
        ppmean(i+1) = nanmean(lvef_class_pp(lvef_class_pp(:,2)==i,1));
    end
    ppmean(1) = NaN;
    % Perform statistical test
    indst = find(lvef_class_pp(:,2) == 0);
    lvef_class_pp(indst,:) = [];
    p_val = kruskalwallis(lvef_class_pp(:,1),lvef_class_pp(:,2),'off');
    hold on
    plot(ppmean,'-og','LineWidth',1.5,'MarkerFaceColor','g','MarkerEdgeColor','k');
    xlabel('LVEF'); ylabel('Pulse Pressure (mmHg)');
    title(sprintf('ICU time: %d hour(s), Total: %d, %d, %d, %d, %d -- %d/%d, p-val = %.3f',j,Npt(1),Npt(2),Npt(3),Npt(4),Npt(5),sum(Npt),sum(nlvef),p_val));
    if DOSAVE
        fileout = sprintf('%slvef_pulsep_mean_tt%d.png',dirout,j);
        print(fig(4),'-dpng',fileout);
    end
    pause;
    hold off
end

%% Analyse SD pulse pressure for different patient categories

for j = test_time
    lvef_class_pp = [];
    IS_SD = 1;
    Npt = zeros(5,1);
    for c = 1: 5
        for i = 1 : length(lvef_class{c})
            ind = find(ptstruct(lvef_class{c}(i)).data{1}.hr(:,1) < j*60);
            if ~isempty(ind)
                vec = ptstruct(lvef_class{c}(i)).data{1}.hr(ind,2);
                if IS_SD
                    vec = std(vec(vec>0));
                else
                    vec = vec(vec>0);
                end
                Npt(c) = 1 + Npt(c);
                lvef_class_pp = [lvef_class_pp; vec (c-1)*ones(size(vec))];
            end
        end
    end
    classLVEF = 0:4;
    ppmean = zeros(1,5);

    fig(5) = figure(102);
    boxplot(lvef_class_pp(:,1),lvef_class_pp(:,2));
    for i = classLVEF
        ppmean(i+1) = nanmean(lvef_class_pp(lvef_class_pp(:,2)==i,1));
    end
    ppmean(1) = NaN;
    % Perform statistical test
    indst = find(lvef_class_pp(:,2) == 0);
    lvef_class_pp(indst,:) = [];
    p_val = kruskalwallis(lvef_class_pp(:,1),lvef_class_pp(:,2),'off');
    hold on
    plot(ppmean,'-og','LineWidth',1.5,'MarkerFaceColor','g','MarkerEdgeColor','k');
    xlabel('LVEF'); ylabel('Pulse Pressure (mmHg)');
    title(sprintf('ICU time: %d hour(s), Total: %d, %d, %d, %d, %d -- %d/%d, p-val = %.3f',j,Npt(1),Npt(2),Npt(3),Npt(4),Npt(5),sum(Npt),sum(nlvef),p_val));
    if DOSAVE
        fileout = sprintf('%slvef_pulsep_sd_tt%d.png',dirout,j);
        print(fig(5),'-dpng',fileout);
    end
    pause;
    hold off
end


%% Analyse time stamps for heart rate and pulse pressure
% Get indices of all patients who have hr and pp data
iAll = [];
for c = 1 : 5
    iAll = [iAll; lvef_class{c}(:)];
end
pt = 1;
ana_timestamps(ptstruct(iAll(pt)).data{1}.hr,ptstruct(iAll(pt)).data{1}.pulsep);

plot_ts(ptstruct(iAll(pt)).data{1}.hr);

X = [];
for i = iAll'
    X = [X min(ptstruct(i).data{1}.pulsep(:,1))];
end

